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Inspired by the parity-time symmetry concept, we show that a judicious spatial modulation of 
gain and loss in epsilon-near-zero metamaterials can induce the propagation of exponentially-bound 
interface modes characterized by zero attenuation. With specific reference to a bi-layer configuration, 
via analytical studies and parameterization of the dispersion equation, we show that this waveguiding 
mechanism can be sustained in the presence of moderate gain/loss levels, and it becomes leaky (i.e., 
radiative) below a gain/loss threshold. Moreover, we explore a possible rod-based metamaterial 
implementation, based on realistic material constituents, which captures the essential features of the 
waveguiding mechanism, in good agreement with our theoretical predictions. Our results may open 
up new possibilities for the design of optical devices and reconfigurable nanophotonics platforms. 
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I. INTRODUCTION 

The possibility to spatially modulate loss and gain 
brings about new dimensionalities in the design of meta- 
materials, which extend far beyond traditional loss- 
compensation schemes. Within this framework, particu¬ 
larly inspiring is the concept of parity-time (VT) symme¬ 
try, originally conceived in quantum physics. 1-3 Against 
the standard assumptions in quantum mechanics, Bender 
and co-workers 1-3 proposed an extended theory where 
the Hermitian property of the Hamiltonian was replaced 
by a weaker symmetry condition on the quantum poten¬ 
tial, V(x) = V*(—x), involving the combined parity (i.e., 
spatial reflection, V) and time-reversal (i.e., complex- 
conjugation, T) operator. They showed that, albeit non- 
Hermitian, such VT -symmetric systems may still exhibit 
entirely real eigenspectra provided that their eigenstates 
are likewise PT-symmetric. However, in view of the an- 
tilinear character of the VT operator, this last condition 
may not hold beyond some non-Hermiticity threshold, 
and the system may undergo a “spontaneous symmetry 
breaking”, i.e., an abrupt phase transition to a complex 
eigenspectrum. 1-3 

In view of the formal analogies between quantum me¬ 
chanics and (paraxial) optics, such concept can be trans¬ 
lated to electromagnetic structures by means of spatial 
modulation of loss and gain, which is becoming tech¬ 
nologically viable. In particular, optical “testbeds” of 
PP-symmetric Hamiltonians have been proposed, 4,5 and 
experimentally characterized in either passive 6 (pseudo- 
PT-symmetric) and actual gain-loss 7 configurations. 
Moreover, a variety of PP-symmetry-inspired exotic ef¬ 


fects have been observed in optical, plasmonic, circuit- 
based, and metamaterial structures, including unidirec¬ 
tional propagation phenomena (invisibility, tunneling, 
negative refraction), coherent perfect absorption, beam 
switching, and absorption-enhanced transmission, with 
very promising potential applications to novel photonic 
devices and components (see Refs. 8-32 for a sparse 
sampling). More recently, potential applications have 
also been proposed in connection with magnetic 33 and 
acoustic 34 structures. It is worth stressing that the po¬ 
tential technical issues that have been recently pointed 
out 35 in connection with the PT-symmetric extension of 
quantum mechanics do not affect these electromagnetic 
and acoustic analogues. 

In this paper, we present a study of PP-synnnetry- 
induced waveguiding in metamaterial slabs. Wave prop¬ 
agation at an interface between two media typically re¬ 
quires one of them to be conducting or with a nega¬ 
tive real part of the permittivity. Here, on the contrary, 
we study the propagation of exponentially-bound modes 
that can be sustained at a gain-loss interface under VT- 
symmetry conditions, without requiring any modulation 
of the real part of permittivity across the interface. This 
intriguing propagation mechanism does not require neg¬ 
ative values of the permittivity (real-part), and it is char¬ 
acterized by a purely real propagation constant (i.e., no 
attenuation). However, in order to achieve substantial lo¬ 
calization along the transverse direction, unfeasibly high 
values of gain are generally required. 11 We therefore sug¬ 
gest to operate in the epsilon near zero (ENZ) regime 36 
(i.e., vanishingly small real part of the permittivities), 
in view of its well-known capabilities to dramatically en- 



2 



FIG. 1. (Color online) Problem schematic. A "PT-symmetric 
bi-layer consisting of two slabs of identical thickness d, and 
relative permittivity distribution as in (1), which can support 
TM-polarized modes exponentially bound and the gain-loss 
interface 2 = 0, and propagating without attenuation along 
the x direction. 


hance the effects of relatively low levels of loss and/or 
gain. 37-39 

Accordingly, the rest of the paper is laid out as follows. 
In Sec. II, we introduce the waveguiding mechanism and 
discuss its attractive features as well as its limitations. 
In Sec. Ill, with specific reference to the ENZ regime, 
we analytically derive the dispersion equation for a VT- 
symmetric bi-layer, and we identify a threshold condition 
on the gain/loss level which separates the bound- and 
leaky-mode regions. In Sec. IV, we explore a possible 
rod-based metamaterial implementation which relies on 
a realistic (semiconductor) gain material. Finally, in Sec. 
V, we provide some concluding remarks and perspectives. 


II. BACKGROUND AND PROBLEM 
STATEMENT 

A. Geometry 


With reference to the schematic in Fig. 1, we start 
considering an isotropic, non-magnetic, piece-wise homo¬ 
geneous bi-layer composed of two slabs of identical thick¬ 
ness d , immersed in vacuum, infinitely extent in the x, y 
plane, and paired along the 2 -direction. Our model is 
hence parameterized by the relative permittivity distri¬ 


bution 


( 1, |s| > d, 

e(z)=l £i, — d < z < 0, (1) 

e*, 0 < z < d, 

where e\ = e' — ie" , with e' > 0, e" > 0. Under the 
assumed time-harmonic [exp(— iujt)\ convention, this im¬ 
plies that the regions —d < z < 0 and 0 < z < d are char¬ 
acterized by gain and loss, respectively, and the structure 
fulfills the necessary condition for VT symmetry, 

£(*) = £*(-*). ( 2 ) 


B. "PT-symmetry-induced surface-waves 

In Ref. 11, it was pointed out that, for transverse- 
magnetic (TM) polarization (i.e., y-directed magnetic 
field), the structure in Fig. 1 may support a 'PT-induced 
surface wave exponentially bound at the gain-loss inter¬ 
face z = 0. This waveguiding mechanism is perhaps 
more easily understood in the half-space limit d —> oo, 
for which the dispersion relationship is simply given by 
(see Appendix A for details) 


k x — ko 


I £i £j 

£l + £i 


= fcn 


(e') 2 + (e"Y 
2e' 


(3) 


with k 0 = w/Co = 27 t/Ao denoting the vacuum wavenum¬ 
ber (and Co and Ao the corresponding wavespeed and 
wavelength, respectively). Accordingly, the field local¬ 
ization in the gain and loss regions is controlled by the 
(complex) transverse wavenumbers 

k z i = \J £i/cq - fc 9 , Im (k z i) < 0, (4) 

and k'* zl , respectively. 

The dispersion relationship in (3) can be interpreted 
as a generalization of the Zenneck-wave 40 and surface- 
plasmon-polariton 41 (SPP) cases, featuring oppositely 
signed imaginary parts of the permittivities. By com¬ 
parison with these two latter cases, the following obser¬ 
vations are in order: 


i) Both media exhibit the same positive value of 
permittivity (real-part), and therefore the mecha¬ 
nism differs substantially from gain-assisted SPP- 
propagation schemes. 42 

ii) The PT-symmetry condition inherently yields a 
real-valued propagation constant k x , i.e., unattenu¬ 
ated propagation along the gain-loss interface. 

in) From the physical viewpoint, such waveguiding 
mechanism is sustained by a transverse (i.e., z- 
directed) component of the power flux from the 
gain- to the loss-region. 
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FIG. 2. (Color online) Geometry as in Fig. 1, but in the 
asymptotic limit d —> oo. Decay length [cf. (5)], scaled by 
the vacuum wavelength, as a function of the gain/loss level s", 
for various representative values of the relative-permittivity 
real part: s' = 10” 4 (squares), s' = 10” 3 (circles), s' = 10” 2 
(up-triangles), s' = 0.1 (down-triangles), s' = 1 (diamonds), 
s' = 10 (stars). 


iv) The branch-cut choice in the gain region [cf. (4)] 
may appear somewhat arbitrary, given that the 
usual radiation condition and decay at infinity can¬ 
not be used as an argument in a gain background. 
Indeed, this a rather controversial issue in the liter¬ 
ature (see, e.g., Refs. 43-48 for a sparse sampling). 
We point out, however, that this choice is irrelevant 
for the bi-layer scenario of actual interest here, and 
it only matters for the half-space configuration. 47 
This latter is, however, an unrealistic limit that 
we consider only in view of the particularly sim¬ 
ple form of the dispersion relationship. Neverthe¬ 
less, for several representative values of s' and e" 
(within and beyond the ENZ regime), we verified 
numerically that the choice in (4) yields results that 
are consistently in agreement with those obtained 
by truncating (along z) the half-space configura¬ 
tion at distances for which the field is sufficiently 
decayed. 

The above waveguiding mechanism looks potentially at¬ 
tractive under many respects. For instance, one may 
envision nanophotonics platforms where channels of gain 
media are suitably embedded in a lossy background, so 
that the waveguiding may be selectively enabled (and 
possibly reconfigured) by optically pumping certain spa¬ 
tial regions. So, effectively we may have “waveguiding on 
demand”, where and when we want it. This may bring 
about new perspectives and degrees of freedom in the de¬ 
sign of optical switches, modulators, and reconfigurable 
photonic networks. 


C. Transverse localization vs. gain/loss level 

Although, in view of (4), the half-space limit always 
features exponential decay (along z) of the fields, one 
intuitively expects the localization to depend critically 
on the gain/loss level (and to vanish in the absence of 
gain and loss). For a more quantitative assessment of 
such localization properties, we show in Fig. 2 the decay 
length 41 

Ld = (5) 

as a function of the imaginary part (absolute value) of the 
permittivity z". for representative values of the real part 
s' spanning several orders of magnitude. As evidenced 
by the log-log scale, for a given wavelength and relative- 
permittivity real-part, the decay length decreases alge¬ 
braically with increasing values of the gain/loss level. In 
particular, localization on subwavelength scales requires 
values of s" that are of the same order or even larger 
than e' . Thus, assuming for instance s' = 10 (compati¬ 
ble with semiconductor materials at optical wavelengths), 
gain/loss levels as high as z" = 3 would be required to 
attain a decay length L d ~ A 0 /4. To give an idea, at the 
telecom wavelength Ao = 1550 nm, this corresponds to a 
gain coefficient 7 = 47rlm (y^el) /Ao ~ 38000 cm” 1 , i.e., 
about an order of magnitude larger than the largest gain 
levels attainable with current technologies. 42,49 ” 51 

What also clearly emerges from Fig. 2 is that de¬ 
creasing e' may allow working with substantially lower 
gain/loss levels. For instance, assuming s' = 10” 4 , de¬ 
cay lengths L d ~ A 0 /4 could be attained with gain/loss 
levels e" ~ 0.009, i.e., gain coefficients (at Ao = 1550nm) 
7 ~ 5000 cm” 1 . 

III. THE ENZ REGIME 

From the above results and observations, it turns out 
that the ENZ regime, 

e' < e" < 1 , ( 6 ) 

seems particularly promising for the waveguiding mech¬ 
anism of interest. While the desired ENZ "PT-symmetric 
characteristics cannot be found in natural materials, we 
show hereafter (see Sec. IV below) that they can be artifi¬ 
cially engineered based on realistic material constituents. 
Before that, however, we study in detail the more realis¬ 
tic bi-layer (i.e., finite d) scenario in Fig. 1 in the ENZ 
regime (6). 

A. Dispersion equation: Bound vs. leaky modes 

It can be shown (see Appendix B for details) that a 
PT-symmetric ENZ metamaterial bi-layer [cf. Fig. 1] 
supports modes propagating along the x direction with a 
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generally complex propagation constant k x which satisfies 
the dispersion equation 

*^o{|n| 2 Re[e?(fc* 1 ) 2 ]-|£i| 2 |^ 1 |}-|fc z i| 2 Re(£iA;* 1 r*) = 0, 

(7) 

where 

k z0 = \Jkl~ k%, (8) 

and 




d/i a 


ri = tan {k z \d). (9) 

In view of the inherent geometrical symmetry, with¬ 
out loss of generality, we focus hereafter on the case 
R e(k x ) > 0 (i.e., propagation along the positive x di¬ 
rection). Among the possible solutions of (7) in the com¬ 
plex k x plane, we are especially interested in bound modes 
characterized by 

Re (k x ) > ko, Im(fc z0 ) > 0, (10) 

i.e., an exponential decay in the exterior vacuum region 
1 2 1 > d. While it is well-known that no such mode can be 
sustained by a low-permittivity slab in the absence of loss 
and gain, it can be shown (see Appendix C for details) 
that this becomes possible for gain/loss levels beyond a 
threshold value 


FIG. 3. (Color online) Geometry as in Fig. 1, for e' = 10 -4 
and d = 0.5Ao. (a) Real- (blue-solid; left axis) and imaginary- 
part (red-dashed; right axis) of the numerically-computed 
[from (7)] propagation constant, as a function of the gain/loss 
level e", illustrating the transition from leaky to bound modes 
occurring at the threshold e" = 0.014 (black-dotted vertical 
line), (b) Gain/loss level threshold [cf. (11)] as a function of 
d/ Ao, for s' = 10 -4 (squares), s' = 10 -3 (circles), s' = 10 -2 
(triangles). 

Figure 3(b) shows the behavior of such threshold as 
a function of the bi-layer electrical thickness, for repre¬ 
sentative values of e'. We observe that the threshold 
depends only mildly on the bi-layer electrical thickness 
and, for sufficiently thick bi-layers (k$d 1, i.e., r 0 ss 1), 
it approaches the asymptotic value 


e'(2 - £')[s'fc 0 d(rg - 1) + 2t 0 ] 
k 0 d(s' - 2)(tq - 1) + 2t 0 


t 0 = tanh (kod ), 


_( n ) 

and it also implies Im(fc x ) = 0 (i.e., no attenuation). 
Below such threshold, leaky modes can instead be found, 
characterized by complex propagation constants 


Re {k x ) < ko, Im (k x ) > 0, Im(fe zO )<0. (12) 


To avoid possible confusion, we stress that the complex 
character of these latter solutions is by no means related 
to the aforementioned spontaneous symmetry breaking 
phenomenon in "PT-symmetric systems, 1-3 as it would 
also arise in the absence of gain and loss. 52 These so¬ 
lutions exhibit exponential decay along the propagation 
direction x, and exponential growth along the transverse 
direction z. Although such character appears clearly un¬ 
physical, they have long been utilized in the antenna com¬ 
munity to effectively model physical resonant radiative 
states in waveguides. 53 

To illustrate the threshold phenomenon, Fig. 3(a) 
shows the numerically-computed propagation constant 
k x as a function of s”, for given values of s' and the 
bi-layer electrical thickness. As it can be observed, for 
increasing values of e" there is a smooth transition from a 
leaky [cf. (12)] to a bound [cf. (10)] mode solution. The 
separation between these two regions occurs at the graz¬ 
ing condition k x = ko , and the corresponding gain/loss 
level is in very good agreement with the analytical esti¬ 
mate of the threshold s" in (11). 


e'L = y/e> (2 - £'), (13) 

which is consistent with enforcing \k x \ > ko in the asymp¬ 
totic dispersion relationship (3). Moreover, as it can be 
expected, the threshold increases with increasing values 
of s', but maintains moderately small values within the 
ENZ regime of interest. We stress that the threshold in 
(11) and its asymptotic limit in (13) are only valid in 
the ENZ limit (6). Therefore, the fact that s in (13) 
vanishes for s' = 2, and it becomes imaginary for s' > 2, 
by no means indicates that the threshold disappears for 
sufficiently thick bi-layers, but rather than the ENZ ap¬ 
proximation is no longer valid in those parameter ranges. 


B. Representative results 

Figures 4 and 5 illustrate some representative results 
for s' = 10 -4 and two feasible gain/loss levels. More 
specifically, for an above-threshold case [s" = 0.02, cf. 
Fig. 3(b)], Fig. 4(a) shows the numerically-computed 
dispersion relationship of a bound mode. As theoretically 
predicted, we observe a purely real propagation constant 
(i.e., no attenuation), which approaches the asymptotic 
prediction [cf. (3)] for d/X 0 > 0.3. To verify the physical 
character of this mode and its actual excitability, Fig. 
4(b) shows a numerically-computed (see Appendix D for 
details) near-field map pertaining to a finite-size (along 
x) structure excited by a magnetic line-source located at 
the gain-loss interface at x = 0. A bound-mode structure 
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FIG. 4. (Color online) (a) As in Fig. 3(a), but as a function of dj Ao (dispersion relationship), for e' = 10~ 4 and e" = 0.02 
(above-threshold case). Also shown (black-dotted horizontal line), as a reference, is the asymptotic limit (3). (b) Numerically- 
computed field magnitude (\H y \) map for a bi-layer with d = 0.5Ao and finite-size (along x ) width of 25Ao (delimited by a 
black-solid rectangle), excited by a magnetic line source located at x = z = 0. Values are sampled so as to avoid the singularity 
at the source, and are normalized with respect to the maximum, (c) Transverse cut (magenta-dashed) at x = 4.17Ao, compared 
with analytical bound-mode prediction [black-solid; cf. (Bl)[ with k x = 1.414fco. 



0 max 


0.30 15 




I.OOu 

0.25 10 



(b) 

0.75- 





CD 0-50 

0.20 5 





nnS" n 


f 


g_0.25 

s > 



' , 

o 0 00 

0.10 -5 




nl 0.25 





"O 

0.05 -10 




£ 0.50 

0.00 -15 




0.75- 

0 

-2 -1 

(J 

1 2 






i.oo-l 


z 

M> 



90 



FIG. 5. (Color online) (a), (b) As in Figs. 4(a) and 4(b), respectively, but for e" = 0.006 (subthreshold case), (c) Numerically- 
computed radiation pattern (with the angle 9 measured with respect to the « axis) compared with leaky-mode-based theoretical 
prediction in (14) for k x = (0.486 + i0.02)fco- 


is clearly visible, with a standing-wave pattern originat¬ 
ing from the structure truncation along the x direction. 
For a more quantitative assessment, Fig. 4(c) shows a 
transverse (z) cut, which clearly exhibits an exponential 
localization, and is in excellent agreement with the the¬ 
oretical prediction (see Appendix B). 


Figures 5(a)-5(c) illustrate the corresponding results 
for a subthreshold gain/loss level [e" = 0.006, cf. Fig. 
3(b)]. More specifically, in the dispersion relationship 
[Fig. 5(a)] we now observe a complex propagation con¬ 
stant, which is indicative of a leaky mode [cf. (12)]. As 
also evident from the near-field map in Fig. 5(b), this 
represents a physical resonant radiative state supported 
by the bi-layer. As a further confirmation, Fig. 5(c) com¬ 
pares the numerically-computed (far-field) radiation pat¬ 


tern with the theoretical leaky-mode-based prediction. 


53 


pr(ff) 

v 


(6) ~ A cos 2 9 ■ 


[fcg sin 2 9 — (/3 2 — a 2 )] + (2 a/3) 2 

(14) 

where k x — /3 + ia is the complex propagation constant 
[cf. (12)], A is a normalization constant, and the angle 9 
is measured with respect to the z axis. A good agreement 
is observed, with the discrepancies attributable to the 
finite-size aperture (along x) of the bi-layer. 


C. Some remarks 

In essence, from a physical viewpoint, the threshold 
phenomenon implies that for low gain/loss levels, the 
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FIG. 6. (Color online) Metamaterial implementation, (a) Unit-cell describing a 2-D array of non-magnetic cylindrical rods of 
radius r c and relative permittivity e c = e' c — is” arranged according to a square lattice with period a. (b) Schematics of the 
effective-medium model: vacuum-coated rod embedded in an effective medium of unknown parameters e e and fi e . The radius 
ro is chosen so that the area of the coated rod is that of the actual square unit-cell in (a), (c) Representative results from 
the synthesis problem in the 3-D parameter space (e", a/Ao, r c /a), assuming e' c = 11.38: each marker represents a candidate 
configuration that satisfies the (asymptotic) condition in (22) for the existence of an unattenuated bound mode. 


transverse power flow from the gain to the loss region 
is not sufficient to sustain a bound mode, and the struc¬ 
ture tends to radiate [at an angle and with a beam-width 
strictly related to the complex propagation constant, 53 
cf. (14)]. This is similar to what is observed in standard 
(lossless, gainless) low-permittivity slabs. 52 By increas¬ 
ing the gain/loss level, the radiation direction progres¬ 
sively departs from the z axis, and becomes grazing at 
the threshold value e" in (11). Beyond this threshold, the 
transverse-power-flow mechanism becomes sufficiently ef¬ 
fective for the structure to sustain a bound mode. 

Incidentally, we found a similar threshold phenomenon 
(with identical parameterization) in a previous study 24 
dealing with the surface-wave-mediated tunneling of im¬ 
pinging waves through the same structure as in Fig. 1. 
This is not surprising, based on reciprocity arguments. 

Another interesting aspect of the above described 
waveguiding mechanism is that the propagation constant 
in the above-threshold (bound-mode) region is inher¬ 
ently real , irrespective of the gain/loss level and elec¬ 
trical thickness. In other words, these bound modes 
are not subject to the spontaneous symmetry break¬ 
ing phenomenon that generally occurs in "PT-symmetric 
systems. 1-3 This is quite unusual, and not observable 
in other waveguiding mechanisms. To give an idea, 
for e' > 1, a PT-symmetric bi-layer could also sup¬ 
port higher-order guided modes, which may be viewed as 
the complex-valued transpositions of the standard guided 
modes supported by a dielectric (lossless, gainless) slab 
waveguide. For such modes, parameters could be tuned 
so as the propagation constant would stay real within cer¬ 
tain ranges. However, by increasing the gain/loss level 
and/or the electrical thickness, spontaneous symmetry 
breaking would eventually occur, and the propagation 


constant would become complex. 


IV. POSSIBLE IMPLEMENTATION 

A typical implementation of ENZ metamaterials at 
optical wavelengths is based on multilayers combining 
thin subwavelength layers of positive- (e.g., dielectric) 
and negative-permittivity (e.g., metals or oxides) ma¬ 
terials. In such implementations, the use of gain has 
been proposed in order to compensate the unavoidable 
loss effects. 54,55 However, since an interface between a 
positive- and negative-permittivity material is naturally 
capable to support a surface wave also in the absence 
of balanced gain and loss, such implementation may not 
allow a clear-cut visualization and interpretation of the 
■P7~-symmetry-induced waveguiding phenomenon of in¬ 
terest here. 

For a more effective illustration of our arguments, we 
therefore take inspiration from all-dielectric implementa¬ 
tions of near-zero-refractive-index metamaterials based 
on periodic arrays of high-permittivity cylindrical rods 
exhibiting Dirac-cone dispersion at the center of the Bril- 
louin zone. 56,57 


A. Effective parameters 

As schematically illustrated in the unit-cell shown in 
Fig. 6(a), we consider a possible implementation con¬ 
sisting of non-magnetic cylindrical rods of radius r c and 
relative permittivity e c = e’ c — ie" c arranged according to 
a square lattice with period a. As in Ref. 56, we model 
such metamaterial by means of the effective-medium the- 
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ory developed in Ref. 58. In essence, as illustrated in 
Fig. 6(b), such model assumes a vacuum-coated cylinder 
of total radius ro embedded in an effective medium of 
unknown parameters e e and /i e . The radius ro is chosen 
so that the area of the coated cylinder is the one of the 
actual square unit-cell, and the effective parameters are 
computed by self-consistency, i.e., by enforcing that the 
total scattering of an electromagnetic wave vanishes. In 
particular, in the limit k e r o <C 1, we obtain 58 the simple 
equations 59 


Ji {kor 0 ) 

fccd’o J[ (fco^o) ; Y{ {k 0 r 0 ) ( D 1 \ 

Y! (fcpro) U[ ( koro ) V1 + A ) ’ 1 ’ 

k 0 r 0 Y{ (koro) 


2 J'q (k 0 r 0 ) 

^ e k 0 r 0 J 0 ( kpro) _ Y 0 (k 0 r 0 ) / Dp \ 
2 Yq (k 0 r 0 ) iJ 0 (kpro) V 1 + D o / 
Ale k 0 r 0 Y 0 (fc 0 r 0 ) 


(16) 


which can readily be solved analytically in closed form. 
In (15) and (16), k c = kp^/e^, J v and Y v are the j/th-order 
Bessel and Neumann functions, 60 respectively, the prime 
denotes differentiation with respect to the argument, and 


_ k< : .7 v ( k c r c ) J v (kpT c ) £ c koJv ( k c v c ) </„ (kpv c ) 

e c kpJ v ( k c r c ) H'} 1] ( k 0 r c ) - k c J' ( k c r c ) ( k 0 r c ) 

(17) 

with denoting the j/th-order Hankel function of the 
first kind. 60 Referring to Ref. 58 for a thorough assess¬ 
ment of the range of applicability of the above model, 
we stress that the underlying approximation does not re¬ 
quire koro, kpr c and k c r c to be small, and thus its validity 
can extend beyond the standard long-wavelength limit. 


B. Model generalizations 


In view of the generally magnetic character of the effec¬ 
tive medium, our VT -symmetric model in Fig. 1 needs to 
be generalized, by assuming also a relative permeability 
distribution 


H(z) 


1, \z\ > d, 

Hi, —d < z < 0, 

fi\, 0 < z < d, 


(18) 


where H\ = // — i\j !', with \j! > 0, n" > 0. Ac¬ 
cordingly, the dispersion relationship of a TM-polarized 
bound mode in the asymptotic (d —> oo) limit can be 
generalized as follows (see Appendix A for details) 


kx — kp\ 


(eiHt - £*Mi) 


= kp |ei| 


4 - tor 

e" H' — e'\x" 
2e"e' ’ 


(19) 


subject to the further condition 

Im j = 0, (20) 

where 

k z i = ijkfeini - k%, Im (Am) < 0. (21) 

For the bi-layer (i.e., finite d) case, the dispersion equa¬ 
tion remains formally identical to (7), but with k z \ de¬ 
fined in (21). In principle, it is also possible to gener¬ 
alize the threshold condition in (11), but the derivation 
is rather cumbersome. Instead, we consider the asymp¬ 
totic limit d/\p » 1 (of direct interest for our subsequent 
studies), for which the existence of a bound mode can be 
established by enforcing in (19) real-valued solutions with 
k x > kp, which yields 

n" fY 2 

-77 < ~T-o-O' ( 22 ) 

e" s' (e") + (e') 

C. Synthesis 

In view of the simple analytical structure of the 
effective-medium model in (15) and (16), and the lim¬ 
ited number of parameters, we found it computationally 
effective to synthesize the metamaterial via a constrained 
parameter search. In what follows, we focus on the syn¬ 
thesis of the gain region, which entails e" > 0; it is easily 
verified from (15) and (16) that the lossy counterpart can 
be obtained by changing the sign of e" c . 

In our synthesis, we fix the real part of the relative per¬ 
mittivity of the rods e' c = 11.38 (compatible with typical 
semiconductor materials), and vary its imaginary part 


0 < e" < 0.35, 

(23) 

the normalized period 


0.1 < a/X 0 < 0.7, 

(24) 

and the normalized cylinder radius 


0 < rja < 0.5. 

(25) 


The above constraints account for the technological fea¬ 
sibility of the required gain level, 42,49-51 the range of va¬ 
lidity of the effective-medium model, 58 and the geomet¬ 
rical consistency of the unit cell, respectively. Figure 
6(c) shows, in the 3-D parameter space (e", a/\p, r c /a), 
a set of possible candidate configurations that satisfy the 
asymptotic condition in (22) for the existence of an unat¬ 
tenuated bound mode. 


D. Results 

As an example, among the possible configurations in 
Fig. 6(c), we consider e" = 0.25, a = 0.465Aq, and r c = 




FIG. 7. (Color online) (a) As in Fig. 4(b), but for a rod-based 
metamaterial implementation with d = 4.65Ao and finite-size 
(along x) width of 27.9Ao. Each half of the bi-layer consists 
of a 10 x 60 square array of cylindrical rods, with period 
a = 0.465Ao, radius r c = 0.375a, and relative permittivity 
s c = 11.38 ± *0.25 (for the gain and loss region, respectively). 
The corresponding effective parameters [cf. (15) and (16)] are 
Eie = 0.002 — *0.107 and /*i e = 0.567 — *0.013 for — d < z < 0 
(gain), and s* e = 0.002 + *0.107 and pi e = 0.567 + *0.013 for 
0 < z < d (loss), (b) Field magnitude {\H y \, normalized with 
respect to the excitation amplitude at a reference plane) at 
the gain-loss interface z = 0 for an infinite (along x) structure 
illuminated by an evanescent plane wave, as a function of the 
k x wavenumber. Also shown as a reference (black-dotted ver¬ 
tical line) is the theoretical bound-mode propagation constant 
[cf. (19)]. 



an unattenuated bound mode with k x = 1.283/co. 

Figure 7(a) shows the numerically-computed field map 
pertaining to the actual rod-based metamaterial struc¬ 
ture excited by a magnetic line source at the gain-loss in¬ 
terface. Also in this case, a bound mode is clearly visible 
and, although the transverse localization is mostly dic¬ 
tated by the microstructure geometry, we can verify that 
the propagation constant is in quantitative good agree¬ 
ment with the theoretical predictions. To this aim, we 
consider an infinite (along x) structure illuminated by an 
evanescent plane wave, and plot in Fig. 7(b) the (normal¬ 
ized) held magnitude at the interface z = 0 as a function 
of the k x wavenumber. We observe that the response is 
strongly peaked around k x = 1.234fc 0 , thereby indicat¬ 
ing a phase-matching with a propagation constant that 
is only ~ 3% different than the theoretical prediction 
above. 

As a further confirmation, we decrease the gain/loss 
level in the rods to e" = 0.05, leaving all other pa¬ 
rameters unchanged. This yields the effective parame¬ 
ters £i e = 0.007 — *0.021 and /* le = 0.567 — *0.003, for 
which the bound-mode condition in (22) is no longer sat¬ 
isfied. Accordingly, numerical solution of the dispersion 
equation (7) [with (21)] now predicts a leaky mode with 
k x = (0.147 + *3.4- 10- 5 )fc 0 . 

Figure 8 shows the results pertaining to the actual rod- 
based structure. In particular, from the field map in Fig. 
8(a) the radiative character of the mode is quite evident. 
Also in this case, looking at the (far-held) radiation pat¬ 
terns in Fig. 8(b) we find a good agreement with the 
theoretical prediction [cf. (14)]. 

Overall, the above results indicate that the rod-based 
metamaterial implementation, based on realistic mate¬ 
rial constituents, reproduces fairly well the waveguid- 
ing mechanism of interest, with good agreement be¬ 
tween numerical simulations and theoretical predictions. 
As previously mentioned, the reliance on material con¬ 
stituents with positive (real-part) permittivity removes 
possible ambiguities on the actual nature of the phe¬ 
nomenon, which can thus be clearly attributed to the 
VT -symmetry. 


FIG. 8. (Color online) (a) As in Fig. 7(a), but for s" = 
±0.05, i.e., Sc = 0.007 — *0.021 and p e = 0.567 — *0.003. 
(b) Numerically-computed radiation pattern (with the angle 
9 measured with respect to the z axis) compared with leaky- 
mode-based theoretical prediction in (14) for k x = (0.147 + 
*3.4 • 10 -5 )fco. 


0.375a, which yields the effective parameters [cf. (15) 
and (16)] £i e = 0.002 — *0.107 and /*i e = 0.567 — *0.013 
for the gain region. Accordingly, the lossy region (£j e = 
0.002 + *0.107 and /u* e = 0.567 + *0.013) can be synthe¬ 
sized by utilizing the same parameters, but e" = —0.25. 
Assuming an idealized VT -symmetric bi-layer with such 
effective parameters and d = 4.65Ao, numerical solu¬ 
tion of the dispersion equation (7) [with (21)] predicts 


V. CONCLUSIONS AND PERSPECTIVES 

To sum up, we have shown that ENZ metamaterial bi¬ 
layers can support "PT-symmetry-induced bound modes 
at the gain-loss interface. These modes propagate with¬ 
out attenuation provided that the gain/loss level exceeds 
a critical threshold, and otherwise exhibit a leaky (ra¬ 
diative) character. Starting from the analytical studies 
and parameterizations, we have designed and simulated 
possible rod-based metamaterial implementations. 

Overall, our results indicate that this intriguing VT- 
symmetry-induced waveguiding mechanism can be ob¬ 
served in the presence of gain/loss levels that are com¬ 
patible with current technological constraints. This may 
set the stage for interesting applications to reconhgurable 
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nanophotonic platforms, as well as novel strategies for the 
design of optical switches and modulators. Besides these 
potential applications, we are currently exploring possi¬ 
ble alternative metamaterial implementations, as well as 
the use of more realistic physical models of gain materi¬ 
als. 


where r/o denotes the vacuum characteristic impedance. 
Finally, by enforcing its continuity at the interface z = 0, 
we obtain 


k z i 

£1 



(A4) 


Appendix A: Details on the asymptotic dispersion 
relationships (3) and (19) 


Assuming the more general (electric and magnetic) sce¬ 
nario of VT -symmetric half-spaces, 


e(z) 


ei, z < 0 , , , _ f n i, z < 0 , 

El, z > 0. ’ \ n\, z > 0, 


.(Al) 

a modal solution exponentially bound at the gain-loss 
interface can be written as 


from which the dispersion relationship in (19) readily fol¬ 
lows by squaring and solving with respect to k x . Note 
that, as a consequence of the squaring, (19) may yield 
spurious solutions which do not satisfy (A4). Hence, the 
additional constraint ( 20 ) [which derives directly from 
(A4)] needs to be enforced. 

The dispersion relationship in (3) immediately follows 
by particularizing (19) to the non-magnetic case /j-| = 1. 


Hy (x, z) = C exp ( ik x x) 


exp (ik z iz ), 
exp {ik zl z ), 


z < 0 , 
z > 0 , 


( A2 ) 

where C denotes a normalization constant, and the con¬ 
tinuity condition at the interface z = 0 is enforced. From 
the relevant Maxwell’s curl equation, we then calculate 
the tangential electric field, 


E x (x, z) 


^ d Jk {x z) 

ikge (z) dz ’ ’ 


(A3) 


Appendix B: Details on the dispersion equation (7) 


For the PT-symmetric bi-layer in Fig. 1, a modal solu¬ 
tion exponentially bound at the gain-loss interface z = 0 
can be expressed as 


Hy (x, z) = exp ( ik x x ) 


Ci exp (- ik z0 z ), 

C 2 exp (ik z iz) + C 3 exp {-ik z \z ), 
C 4 exp {ik*\z) + C 5 exp {—ik* zl z ), 
C 6 exp (ik z0 z ), 


z < -d, 

—d < z < 0, 
0 < z < d, 
z > d, 


(Bl) 


with k z i and k z 0 given in (4) and ( 8 ), respectively, and 
the unknown expansion coefficients Cj, j = 1, ...,6 to be 
calculated by enforcing the continuity of the magnetic 
[(Bl)] and electric [cf. (A3) with (1)] tangential fields 
at the three interfaces z = 0 and z = ±d. This yields 
a 6 x 6 homogeneous linear system of equations, whose 
nontrivial solutions can be found by zeroing the system- 
matrix determinant, viz., 


det = e\k1iTi (k zl — ieXk z oT 3 ) 

+ E\k z0 k* zl Ti ( £*ik z0 - ik* z \Ti) 


+ S\k z \ [ 2 ie\k z0 k* zl + (^) 2 k 2 z0 rX + (, k* zl f r* 

= ik z0 |ki | 2 Re e\ (k* zl ) 2 - |£i | 2 |fc z i| 2 | 

- \k z i\ 2 Re {s\k* zl Tl) - Re (\ei\ 2 £ik 2 z0 k* zl T^J (B2) 


neglecting [in view of the assumed ENZ regime, cf. ( 6 )] 
the third-order term in £\. 


Appendix C: Details on the leaky-to-bound mode 
transition 

We now prove that, for gain/loss levels beyond the 
threshold e" in (11) the VT -symmetric ENZ bi-layer in 
Fig. 1 supports a bound mode [cf. (10)] with real prop¬ 
agation constant (i.e., no attenuation). To this aim, it is 
expedient to rewrite the dispersion equation (A4) as 

F (k x ) = 0, (Cl) 


with 


where the last equality follows from simplifications ex¬ 
ploiting the VT -symmetric character. The dispersion re¬ 
lationship in (A4) readily follows by zeroing (B2) and 


F(k x ) = -a 0 {|ri| 2 Re[£?(fc; i ) 2 ] - M 2 |fc 2 i|} 

- IM'Re^rTi*), (C2) 
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and 


a 0 = —ik z o = \Jkl - k %. (C3) 

In such a way, the real character of the dispersion equa¬ 
tion in the parameter range of interest k x > ko is empha¬ 
sized, and a simple bracketing strategy can be exploited 
to prove the existence of real-valued roots. 

First we consider the asymptotic limit k x ko, for 
which we straightforwardly obtain from (C3) and (4) 

a °L>fc 0 ~ ^ kzl \k x »k 0 lkx ( C4 ) 
and hence, from (9), 


T\ 


k x k 0 


(C5) 


By substituting (C4) and (C5) in (C2), we then obtain 


F(k x ) 


k x k 0 


k l £ ' + R- e ( £ l) + kil 2 

e'kl > 0, 


(C6) 


where the last approximate equality stems from neglect¬ 
ing [in view of the assumed ENZ regime, cf. (6)] second- 
order terms in £\. We have thus shown that the left-hand- 
side of the dispersion equation (Cl) is always positive in 
the asymptotic limit k x ^> ko- 

Next, we consider the grazing condition k x = ko, for 
which (C3), (4) and (9) yield 


&o 



= 0, 



and 


k 0 V £ i - 1 



Tl 


k x — k 0 


,£\kod 
—tTO — l - 


•0 ’ 


(C8) 


respectively, with the approximate equality stemming 
from first-order McLaurin expansions in £i. Substitution 
of (C7) and (C8) in (C2) finally yields 


F(ko) 






+ \£x\ 2 k 0 d(£' - 2) (tq - l) 


(C9) 


Recalling the asymptotic behavior in (C6), we can con¬ 
clude that if F (ko) < 0, the dispersion equation in (C2) 
must admit a real-valued solution k x > ko, which corre¬ 
sponds to a bound mode [cf. (10)]. By solving (C9) with 
respect to the gain/loss level e", this condition can be 
parameterized as 


£ > £ f 


(CIO) 


with the threshold e" given in (11). Moreover, since it can 
be numerically verified that, within the parameter range 
of interest, F(k x ) is a monotonic function, the above con¬ 
dition turns out to be not only sufficient, but also neces¬ 
sary. 

For subthreshold gain/loss levels, complex-valued so¬ 
lutions are found instead, which generally exhibit the 
leaky-mode character in (12). 

Appendix D: Details on the numerical simulations 


All the numerical simulations in our study are carried 
out by means of the finite-element-based commercial soft¬ 
ware package COMSOL Multiphysics. 61 In particular, we 
utilize the RF module and the frequency-domain solver. 

For the finite-size configurations in Figs. 4(b), 5(b), 
7(a) and 8(a), we utilize a magnetic line-current excita¬ 
tion located at the center of the structure (x = z = 0), 
perfectly-matched-layer terminations for the computa¬ 
tional domain, and a triangular mesh with adaptive ele¬ 
ment size. This results in a number of elements on the 
order of 2.8 • 10 5 and 1.3 • 10 6 for the idealized [cf. Figs. 
4(b) and 5(b)] and rod-based [cf. 7(a) and 8(a)] config¬ 
urations, respectively. The (far-held) radiation patterns 
[cf. Figs. 5(c) and 8(b)] are straightforwardly obtained 
by utilizing the post-processing tools in the RF module. 61 

The results in Fig. 7(b) refer instead to an infinite 
(along x ) structure, simulated by means of a unit-cell 
consisting of a single row of rods with phase-shift bound¬ 
ary conditions, and excited via a wave-port 61 by an 
evanescent plane-wave. 
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